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Wavelets and Fast Numerical Algorithms 
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Numerical algorithms using wavelet bases are similar to other transform methods 
in that vectors and operators are expanded into a basis and the computations take place 
in the new system of coordinates. As in all transform methods, such approach seeks an 
advantage in that the computation is faster in the new system of coordinates than in the 
original domain. However, due to the recursive definition of wavelets, their controllable 
localization in both space and wave number (time and frequency) domains, and the 
vanishing moments property, wavelet based algorithms exhibit a number of new and 
important properties. 

In the usual transform methods, the functions of the basis (e.g. exponentials, 
Chebyshev polynomials, etc.) are chosen to be eigenfunctions of some differential opera- 
tor (e.g. solutions of the Sturm-Liouville problem). The choice of the differential operator 
and, hence, of the basis functions, is dictated by the availability of fast algorithms for 
expanding an arbitrary function into the basis. Unfortunately, classes of operators which 
have a sparse representation in such bases are very narrow. 

Wavelets, on the other hand, are not solutions of a differential equation. These 
functions are defined recursively and are generated via an iterative algorithm. They are 
translations and dilations of a single function Instead of diagonalizing some differential 
operator, representations in the wavelet bases reduce a wide class of operators to a sparse 
form. Here the orthogonality of wavelets to the low degree polynomials (the vanishing 
moments property) plays a crucial role 0. 

The orthonormal bases of wavelets were fisrt constructed by Stromberg and 
then by Meyer p5| . Later, the notion of the Multiresolution Analysis was introduced by 
Meyer and Mallat P3|. Orthonormal bases of compactly supported wavelets were 



constructed by Daubechies [T^ . There are many new constructions of orthonormal bases 
with a controllable localization in the time-frequency domain, notably "wavelet-packet" 
bases in [0 and |15|, local trigonometric bases in |[T^ and ||2^, wavelet bases on the 
interval in [|l2l and Very important connection exists between the wavelets 



and the technique of subband coding in signal processing. In fact, the discrete wavelet 



^It is possible to construct bases with translations and dilations of several functions, see e.g. 
■^This property and the fact that the basis is orthonormal, distinguish the wavelet bases from the 
hierarchical bases. 
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transform is accomplished by the pair of the so-called quadrature mirror filters. The 
exact quadrature mirror filters (QMFs) were introduced by Smith and Barnwell P^ . 

Wavelets have some of their historical roots in Littlewood-Paley and Calderon- 



Zygmund theories (see e.g. ||2^) which has been powerful tools in analysis of linear 



and non-linear operators. In Numerical Analysis some of the ingredients of Calderon- 
Zygmund theory appear in the Fast Multipole Method (FMM) for computing potential 
interactions |3^, [|19|, FMM was designed for computing potential interactions 



between N particles in 0(— logeA^) operations (instead of 0{N ) operations). The re- 
duction of the complexity in FMM is achieved by approximating the far field effect of 
a cloud of charges located in a box by the effect of a single multipole at the center of 
the box. All boxes are then organized in a dyadic hierarchy enabling an efficient 0{N) 
algorithm. 

Fast wavelet-based algorithms of H| provide a systematic generalization of the 



FMM and its descendents (e.g. |29], B], to all Calderon-Zygmund and pseudo- 



differential operators. The subdivision of the space and its organization in a dyadic 
hierarchy are a consequence of the multiresolution properties of the wavelet bases, while 
the vanishing moments of the basis functions make them useful tools for approximation. 

A novel aspect of representing operators in the wavelet bases is the so-called non- 
standard form P]. The remarkable feature of the non-standard form is the uncoupling of 
the interactions between the scales. The non-standard form leads to an order algorithm 
for evaluating operators on functions. It is also quite remarkable that the error estimates 
for the non-standard form lead to a proof of the selebrated "T(l)" theorem of David and 
Journe (see P]). The non-standard forms of many basic operators, such as derivatives, 
fractional derivatives, the Hilbert and Riesz transforms, may be computed explicitly [Q. 
A straightforward realization, or the standard form, by contrast, contains matrix entries 
reflecting "interactions" between all pairs of scales. The standard form yields, in general, 
only an order A^log(A^) algorithm for evaluating operators on functions. 

The representation of wide classes of operators in wavelet bases may be viewed as 
a method for their "compression", i.e., conversion to a sparse form. For these operators 
sparse representations lead to fast algorithms for matrix multiplications. Since the per- 
formance of many algorithms requiring multiplication of dense matrices has been limited 
by 0{N^) operations, these fast algorithms address a critical numerical issue. 

Among the algorithms requiring multiplication of matrices is an iterative algo- 



rithm for constructing the generalized inverse |31|, the scaling and squaring method for 



computing the exponential of an operator, and similar algorithms for sine and cosine of 
an operator, to mention a few. By replacing the ordinary matrix multiplication in these 
algorithms by the fast multiplication in the wavelet bases, the number of operations 
is reduced to, essentially, an order operations. For example, if both, the operator 
and its generalized inverse, admit sparse representations in the wavelet basis, then the 
iterative algorithm pT| for computing the generalized inverse requires only 0{NlogK) 
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operations, where k is the condition number of the matrix. Various numerical examples 
and applications may be found in [0, |T|] and p. 

Solving the two-point boundary value problem for the elliptic differential opera- 
tors in the wavelet "system of coordinates" allows us to construct the Green's function 
(the inverse operator) in 0{N) operation. We note that the ordinary matrix represen- 
tation of the Green's function requires 0{N'^) significant entries but the representation 
of the Green's function in the wavelet bases requires (for a given accuracy) only 0{N) 
entries. The main tool in constructing the Green's function numerically is the diagonal 
preconditioner available for the periodized differential operators in the wavelet bases 0, 
§ (see also H). 

Unfortunately, the format of one lecture does not allow us to cover all the de- 
velopments or mention all the results available today. Instead, we will review several 
features of the new numerical methodology based on the wavelet representations. Start- 
ing from the notion of multiresolution analysis, we will consider the non-standard form 
(which achieves decoupling among the scales) and the associated fast numerical algo- 
rithms. Examples of non-standard forms of several basic operators (e.g. derivatives) will 
be computed explicitly. 
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I Mult iresolut ion Analysis and Wavelets. 



We briefly outline here the properties of compactly supported wavelets and refer for 
details to [|1^, and [Q. Let us start with the notion of the multiresolution analysis 



26| , which captures the essential features of all multiresolution approaches developed 



so far. 

Definition I.l Multiresolution analysis is a decomposition of the Hilbert space L^(R'^), 
d > 1, into a chain of closed subspaces 

...C V2C ViC VoCV_iCV_2C... (1.1) 

such that 

^- Cljez^j = {0} and Ujez^j ^-5 dense in L^(R'^) 

2. For any f G L^(R^) and any j G Z, f{x) G Vj if and only if f{2x) G Vj_i 

3. For any f G L^(R*^) and any k G Z^, f{x) G Vq if and only if f{x — A;) G Vq 

4. There exists a scaling function ip G Vq such that {ip{x — A;)},tezd is a Riesz basis of 
Vo. 

In this lecture we use only orthonormal bases, so that we replace Condition 4 by 

4'. There exists a scaling function if G Vq such that {^{x — k)}k^2^d is an orthonormal 
basis o/Vq. 

Let us define the subspaces as an orthogonal complement of in Vj_i, 

V,_i = V,©W„ (1.2) 
and represent the space L^(R*^) as a direct sum 

L2(Rd) = 0w,. (1.3) 



Selecting the coarsest scale n, we may replace the chain of the subspaces (|1.1|) by 

V„ C . . . C V2 C Vi C Vo C V_i C V_2 C . . . , (1.4) 

and obtain 

L2(R'1)=V„0W,-. (1.5) 



If there is a finite number of scales then, without loss of generality, we set j = to be 
the finest scale and consider 



V„ C . . . C V2 C Vi C V 
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:i.6) 



instead of (|1.4| ). In numerical realizations the subspace Vq is finite dimensional. 

The function Lp is the so-called scaling function and, with its help, we may define 
the function t/^, the wavelet, such that the set of functions A;)}fegz is an orthonormal 

basis of Wq, 



An example of the multiresolution analysis satisfying Definition |]T| with Condi- 
tion 4' is the chain of subspaces generated by the Haar basis [^. The scaling function in 
this case is the characteristic function of the interval (0, 1). The Haar function is defined 

as 

r 1 for < a; < 1/2 
h{x) = \ -1 for 1/2 < a; < 1 , (1.7) 
[ elsewhere. 

and the Haar basis is formed by functions hj^^{x) = 2^^^'^h{2~^ x — k) j, k E Zi. 

The wavelet bases (with a smooth scaling function of Condition 4') generalizing 
the Haar functions were first constructed by Stromberg and then Meyer p5[. The 

"si and 



notion of the Multiresolution Analysis was introduced by Meyer |^ and Mallat 

p5[] and, of course, of |^ 



it is more recent than the constructions of [g3], and, of course, of |gO[- Compactly 
supported wavelets with vanishing moments were constructed by I. Daubechies [|l^ and 
we will use them in this lecture. However, most of the results that we discuss do not 
depend on this particular choice of the wavelet bases. 

The vanishing moments property simply means that the basis functions are chosen 
to be orthogonal to the low degree polynomials, namely, if the set of functions {iIj{x — 
k)}k£Z is an orthonormal basis of Wq, then 



+00 



ip{x)x"^dx = 0, m = 0, 



M- 1. 



For the Haar function in (|1.7|) M = 1 and it is trivially orthogonal to constants. 



There are two immediate consequences of Definition |]1| with Condition 4'. First, 
the function ip may be expressed as a linear combination of the basis functions of V_i. 
Since the functions {ipj^k{x) = 2~^/'^Lp{2~^x — k)}k&7. form an orthonormal basis of Vj, 
we have 



L-l 



(p{x) = a/2 ^ hk(p{2x — k). 



:i.9) 



k=0 



In general, the sum in ( |1.9| ) does not have to be finite and, by choosing the finite sum in 
L.9|), we are selecting the compactly supported wavelets. We may rewrite (|1.9|) as 



m = m,{i/2)^{i/2), 



;i.io) 
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where 

^(0 = -^/ ^ix)e'^^dx, (1.11) 



oo 



/27r J- 

and the 27r-periodic function mo is defined as 

L-l 



^o(0 = 2-'/' E ^^e'^^- (1-12) 

fc=0 

Second, the orthogonahty of {(p{x — k)}k£z imphes that 

/+00 f + OO 

if{x-kMx)dx= / \ipiO\^e-'^U(, (1.13) 
oo J —oo 



and, therefore, 

r2iT 



4o= /^El^(e + 2vr/)pe-^'=«tie, (1-14) 
and 

Y^m^ + 2nl)\' = l. (1.15) 



Using ( |1.1CI| ), we obtain 

^|mo(e/2 + vr/)n^(e/2 + vr/)|2 = l, (1.16) 
and, by taking the sum in (|1.16| ) separately over odd and even indeces, we have 



E |mo(e/2+27r/)n^(e/2+27r/)|2+^ \mo{^/2+2nl+n)\^\0{^/2+27rl+7r)\^ = 1. (1.17) 
lez /gz 

Using the 27r-periodicity of the function mo and (|1.15|) , we obtain (after replacing ,^/2 by 
^) a necessary condition 

|rrio(Ol' + l"^o(e + vr)|' = l, (1.18) 
for the coefficients hj. in ( |1.12| ). On denoting 



mi(0 = e-'«mo(e + vr), (1.19) 

and defining the function ip, 

iP{x) = V2Y,gkif{2x-k), (1.20) 

k 

where 

gk = i-lfhL-k-i, A; = 0,...,L-1, (1.21) 



or, the Fourier transform of ip, 



m = mi(e/2)^(e/2), (1.22) 

where 

mi (0 = 2-1/2 '"^"'^.e'^^C, (1.23) 

fc=0 



it is not difficult to show (see e.g., |2S|, |T^, [0), that on each fixed scale j G Z, the 
wavelets {i/jj^ki^) = 2^^^'^'ip(2^^x — k)}k,^z form an orthonormal basis of Wj. 



Equation (|1.18|) defines a pair of the quadrature mirror filters (QMFs) H and G, 



where H = {hk}tzl^'^ and G = {gk}lzl^'\ The exact QMF filters were first introduced 



by Smith and Barnwell ||32| for subband coding. 

We will not go into the details of considering necessary and sufficient conditions 
for the quadrature mirror filters H and G to generate the wavelet basis and refer to 
for the details. The coefficients of the quadrature mirror filters H and G are computed by 
solving a set of algebraic equations (see e.g. [|l3)- The number L of the filter coefficients 



in ( |1.12| ) and ( |1.23D is related to the number of vanishing moments M , and L = 2M 



for the wavelets constructed in |T^. If additional conditions are imposed (see [§] for an 



example), then the relation might be different, but L is always even. 

We observe that once the filter H has been chosen, it completely determines 
the functions and ip and therefore, the multiresolution analysis. Moreover, in properly 
constructed algorithms, the values of the functions ip and are (almost) never computed. 
Due to the recursive definition of the wavelet bases, all the manipulations are performed 
with the quadrature mirror filters H and G, even if they involve quantities associated 
with (p and tp. 

As an example, let us compute the moments of the scaling function (j). The 
expressions for the moments. 



MZ = J x"" ip{x) dx, m = 0,...,M-l, (1.24) 
may be found in terms of the filter coefficients {hk}lZi using 

oo 



where mQ{^) is given in (|1.12|) . 



The moments A^™ are obtained (within the desired accuracy) by recursively gen- 
erating a sequence of vectors, {Ai^}m=o'~^ for r = 1, 2, . . . , 



' m 



^T+i = E 2-''-M';-^Mi, (1.26) 
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starting with 

k=L-l 

MT = 2-"'-^ hkk"", m = 0,...,M -1. (1.27) 

Each vector {Ai^}m=o^~^ represents M moments of the product in ( |1.25| ) with r terms, 
and the iteration converges rapidly. Notice, that we never computed the values of the 
function if itself. 
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II The non-standard form 

The wavelet bases in L^(R'^), d > 2, may be constructed as a tensor product of the 
one-dimensional bases. Considering d = 2 and using the Haar basis as an example, 
we note that the supports of the basis functions are rectangles of various dyadic sizes. 
Representing operators in such bases leads to the standard form which we will discuss in 
the next Section. 

Alternatively, wavelet bases in L^(R'^), d > 2 may be constructed using the 
scaling function in addition to the wavelets. Such construction is specific to wavelet 
bases. Considering c? = 2 as an example, we note that the triplet of functions 

{i^j,ki^)i^j,k'iy), ^pj,kix)ipj^k'iy), (Pj,kix)^j,k'iy)}, (2.1) 

where j,k,k' G Z, forms the basis of L^(R^). We note that the basis functions have 
square supports. Representing operators in these bases leads to the non-standard form 

i- 

Let us introduce the non-standard form in the context of the Multiresolution 
Analysis, independently of the specific choice of the wavelet basis. Let T be an operator 

T : V{R) V{R), (2.2) 

with the kernel K{x,y). We define projection operators on the subspace Vj, j G Z, 

Pj:V{R)^Vj, (2.3) 

as follows 

iP,f)ix) = Y.{f,cp^,,)cp,4x). (2.4) 

k 

Expanding T in a "telescopic" series, we obtain 

T = J2(QiTQj + QjTPj + PjTQj), (2.5) 

where 

Q, = - P, (2.6) 

is the projection operator on the subspace W^. If there is the coarsest scale n, then 
instead of (|275| ) we have 

n 

T= (QjTQj + QjTP, + P,TQ,) + PnTP^, (2.7) 

j=-oo 

and if the scale j = is the finest scale, then 

n 

To = T^iQjTQj + QjTPj + PjTQj) + Pr^TPn, (2.8) 
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where T ~ Tq = PqTPq is a discretization of the operator T on the finest scale. 

The non-standard form is a representation (see [Q) of the operator T as a chain 
of triplets 

T = {A„i?„r,},ez (2.9) 
acting on the subspaces Vj and Wj, 

Aj ■ W, ^ W^-, (2.10) 

Bj-.Yj^Wj, (2.11) 

r, : W, ^ V,. (2.12) 

The operators {Aj, Bj, Tj}j^z are defined as Aj = QjTQj, Bj = QjTPj and Tj = PjTQj. 
These operators admit a recursive definition via the relation 

where operators Tj = PjTPj, 

T,- : V, ^ V,. (2.14) 
If there is a coarsest scale n, then 

T = {{Aj, Bj, ^j}jeZ:j<n, Tn}, (2.15) 

where T„ = PnTPn. If the number of scales is finite, then j = 1, 2, . . . , n in ( p.l5|) and 
the operators are organized as blocks of the matrix (see Figures 1 and 2). Let us make 
the following observations: 

1) . The map ( |2.10| ) implies that the operator Aj describes the interaction on the scale j 
only, since the subspace is an element of the direct sum in ( |1.5|) . 

2) . The operators Bj, Tj in ( |2.11| ) and ( ^.12 ) describe the interaction between scale j 



and all coarser scales. Indeed, the subspace Vj contains all the subspaces V^v with j' > j 
(see (O))- 

3). The operator Tj is an "averaged" version of the operator T^-i. 

The operators Aj, Bj and Tj are represented by the matrices , (3^ and 7-', 



<^k,k' = J J K{x,y)iljj^kix)^j,k'iy) dxdy, (2.16) 
PLk' = K{x,y)ipj^k{x)(fij,k'{y) dxdy, (2.17) 



and 

7lfc' = J J K{x,y)(pj^kix)ipj^k'iy) dxdy. (2.18) 
The operator Tj is represented by the matrix s\ 

sik'= K{x,y)ipj^k{x)^j^k'{y) dxdy. (2.19) 
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Figure 1: Organization of the non-standard form of a matrix. The submatrices Aj, Bj, 
and Tj, j — 1, 2, 3, and Ta are the only non-zero submatrices. 
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12 




Figure 3: The non-standard form of the same matrix as in Figure |^ in the basis of the 
wavelets on the interval 0. The vertical and horizontal bands (which are present in 
Figure ^ due to periodization) do not appear in this representation 
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Ill The standard form 



The standard form is the representation of an operator in the tensor product basis. 
Instead of introducing the standard form in this manner, we emphasize the connection 
with the non-standard form. The standard form is obtained by representing 

V, = 0W,., (3.1) 

3'>3 

and considering for each scale j the operators {-Bj , F^' }j'>j, 

5j':W,v ^W,-, (3.2) 

rf :W,^W,,. (3.3) 
If there is the coarsest scale n, then instead of (|3.1|) we have 

j'=n 

V, =V„ W,.. (3.4) 

In this case, the operators {-Bj , F^ } for j' = j + 1, . . . , tt, are as in (^]2|) and (|3.3| ) and, 
in addition, for each scale j there are operators {B^^^^} and {F^"*"^}, 

: V„ ^ W„ (3.5) 

r^"+i : ^ V„. (3.6) 

(In this notation, FJ^+^ = F„ and -B^"*"^ = -Bn). If there are finitely many scales and Vq 
is finite dimensional, then the standard form is a representation of Tq = PqTPq as 

^0 = {Aj, {B^j {F^- };]v=J+i, Bj~^^, Tn}j=i,...,n- (3.7) 

The operators ( |3.7| ) are organized as blocks of the matrix (see Figures 3 and 4). 

If the operator T is a Calderon-Zygmund or a pseudo-differential operator then, 
for a fixed accuracy, all the operators in ( |3.7| ) (except T„) are banded. As a result, the 
standard form has several "finger" bands which correspond to the interaction between 
different scales. For a large class of operators (pseudo-differential, for example), the 
interaction between different scales characterized by the size of the coefficients of "finger" 
bands, decays as the distance j' — j between the scales increases. Therefore, if the scales 
j and j' are well separated, then for a given accuracy, the operators Bj , F^- can be 
neglected. 

There are two ways of computing the standard form of a matrix. First consists in 
applying the one-dimensional transform to each column (row) of the matrix and, then, 
to each row (column) of the result. Alternatively, one can compute the non-standard 
form and then apply the one-dimensional transform to each row of all operators B^ and 
each column of all operators F^ . We refer to for details. 
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Figure 4: Organization of a matrix in the standard form 




Figure 5: An example of a matrix in the standard form (see Example 1) 
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IV Compression of operators 

If the operator T is a Calderon-Zygmund or a pseudo-differential operator then, by using 
the wavelet basis with M vanishing moments, we force operators {Aj, Bj^ Fjjjgz to decay 
roughly as l/rf*^"'"^, where (i is a distance from the diagonal. For example, let the kernel 
satisfy the conditions 

\K{x,y)\ < (4.1) 
F - y\ 

\d^K{x,y)\ + \dl'K{x,y)\ < , j^' (4.2) 

\x y\ 

for some M > 1. Then by choosing the wavelet basis with M vanishing moments, the 
coefficients ali, [3li,'~^li of the non-standard form (see ( p.l6| ) - ( p.l8|) ) satisfy the estimate 



H^ + \PU + H^<Y^0J^^ (4-3) 

for all 

\i-l\>2M. (4.4) 



If, in addition to (|A]), (p|), 

I / K{x,y) dxdy \< C\I\ (4.5) 

JIy.1 



for all dyadic intervals / (this is the "weak cancellation condition", see [^), then ( [4.3|) 
holds for all i, /. 

If T is a pseudo-differential operator with symbol cT(a;, ^) defined by the formula 
T(/)(x) = a{x,D)f = J e^^« a(x,0/(0 = J K{x,y)f{y) dy, (4.6) 

where K is the distributional kernel of T, then assuming that the symbols a of T and 
a* of T* satisfy the standard conditions 

I 5f a(x,0 |< C^A^+ I e 1)'-"+^ (4.7) 
I 9f a*{x,0 \< Q/3(l+ I e l)'^""'^ (4.8) 



we have the inequality 



for all integer 
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Suppose now that we approximate the operator Tq by the operator T(f obtained 
from To by setting to zero all coefficients of matrices , (3^ and 7-' outside of bands of 
width B > 2M around their diagonals. We obtain 

||To^-To||<-^log2iV, (4.10) 

where C is a constant determined by the kernel K and log2 N is the number of scales 
in the representation. In most numerical applications, the accuracy e of calculations is 
fixed, and the parameters of the algorithm (in our case, the band width B and order M) 
have to be chosen in such a manner that the desired precision of calculations is achieved. 
If M is fixed, then we choose B so that 

i?> (^-log^iVj . (4.11) 

In other words. To has been approximated to precision e with its truncated version, which 
can be applied to arbitrary vectors for a cost proportional to N ({C/e) log2 N)^^^\ which 
for all practical purposes does not differ from A^. 

A more detailed investigation permits the estimate ([4.10|) to be replaced with 
the estimate 

l|ro^-ro||<-^, (4.12) 

making the application of the operator To to an arbitrary vector with arbitrary fixed 
accuracy into a procedure of order A^. Obtaining this uniform estimate leads to a proof 
of 

Theorem (G. David, J.L. Journe) Suppose that the operator 

T(/)= f K{x,y) f{y) dy (4.13) 



satisfies the conditions ( [4.1| ), ( [4.2| ), ( |4.5| ). Then a necessary and sufficient condition for 
T to be bounded on is that 

/3(x) = T(l)(x), (4.14) 
7(2/) = T*(l)(y) (4.15) 
belong to dyadic B.M.O., i.e. satisfy condition 

sup f \(3{x) - mj{p)\^dx < C, (4.16) 

where J is a dyadic interval and 

mj{(3) = 7^ f (3{x)dx. (4.17) 



\J\ Jj 
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Again we refer to for details. 

The compression of operators results in fast algorithms for evaluation of operators 
on functions. We present here one example and refer to P] for additional examples. 

Example 1. 

In this example, we consider the matrix 

Aij = 

and convert it to the non-standard form using wavelets with six vanishing moments. 
Setting to zero all entries whose absolute values are smaller than 10"^, we obtain the 
non-standard form where the non-zero elements are shown in black in Figure |^. The 
results of experiments in applying this sparse matrix to a vector are tabulated in Table 1. 
The standard form of the operator A with = 256 is depicted in Figure ^. 

Column 1 of Table 1 contains the number N indicating the size of N x N matrix 
Aij, columns 2, 3 contain CPU times Tg, required by the standard order 0(A^^) 
and the fast 0{N) schemes to multiply a vector by the matrix, and column 4 contains 
the CPU Td time used to produce the non-standard form of the operator. Columns 5, 6 
contain the L2 and Loo errors of the direct calculation, and columns 7, 8 contain the same 
information for the result obtained by computing in the wavelet system of coordinates. 
Finally, the last column contains the compression coefficients Ccomp, defined by the ratio 
of A^^ to the number of non-zero elements in the non-standard form of of the matrix. 



Input 
Size 


Time 


Error of Single Precision 
Multiplication 


Error of FWT 
Multiplication 


Compression 
Coefficient 


N 


Ts 


T 


Td 


L2 - norm 


Loo - norm 


L2 - norm 


Loo - norm 


Ccomp 


64 


0.12 


0.16 


7.76 


1.26- 10-^ 


3.65 • 10-^ 


8.89 • 10-8 


1.72- lO-'^ 


1.39 


128 


0.48 


0.38 


32.62 


2.17- 10"^ 


8.64- lO"'^ 


1.12 • 10"^ 


9.94- lO"'^ 


2.22 


256 


1.92 


0.80 


96.44 


2.81 • lO-'^ 


1.12 • IQ-^ 


1.25 • 10-^ 


5.30- lO-'^ 


3.93 


512 


7.68 


1.80 


252.72 


4.21 • lO-'^ 


1.75- 10"6 


1.23-10-^ 


5.16- lO-'^ 


7.33 


1024 


30.72 


3.72 


605.74 


6.64- 10-^ 


3.90-10-6 


1.36-10-^ 


5.04- lO-'^ 


14.09 



Table 1: Numerical results for Example 1 
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and 
where 



V The operator d/dx in wavelet bases. 

For a number of operators (e.g., differential operators, fractional derivatives, Hilbert and 
Riesz transforms) we may compute the non-standard form in the wavelet bases by solving 
a small system of linear algebraic equations @]. As an example, we construct the non- 
standard form of the operator d/dx. The matrix elements a^^, and of Aj, Bj, and 
Tj, where i, /, j G Z for the operator d/dx are easily computed as 

/•oo 

5.1 
5.2 

5.3 

5.4 
5.5 

5.6 

5.7 

5.8 

5.9 



and 



all = 2"-' / ip{2-^x - i) iIj'{2~^x - I) 2-^dx = 
J —00 

/oo 
i;{2-^x - i) ip'{2-^x - /) 2-^dx = 2-^'A_;, 
-00 

= 2-^ / ^{2-^x - i) i)'{2-^x - I) 2-^dx = 2-^7i_i, 



30 

ai = I iplx — I) —il){x) dx, 
-00 dx 

r+00 d 

/3i = ij{x - I) —(p{x)dx, 
3 dx 



r+00 d 

'yi = I (p{x - I) — V'(x) dx. 
3 dx 



Moreover, using ( |1.9|) and ( [1.20| ) we have 

L-l L-1 



and 



where 



= 2 ^ Qk gy r2i+k-k' , 

k=0 k'=0 
L-l L-l 

A = 2 ^ ^ Qkhki r 2i+k-k', 

k=0 k'=0 
L-l L-l 

7j = 2 X! X! Ok' r2i+k-k', 

k=0 k'=0 

'■+°° , d 



n 



/ + 00 H 
f{x-l)—(f{x)dx, leZ. (5.10 
-00 dx 



Therefore, the representation of d/dx is completely determined by the coefficients ri in 
( p.lU| ) or in other words, by the representation of d/dx on the subspace Vq. Rewriting 
( p.lOp in terms of (p{^) (see ( |1.11|) ), we obtain 

r+00 

ri= \mm>-''^d^- (5.11) 
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Thus, the coefficients r; depend only on the autocorrelation function of the scaling func- 
tion yj, rather that the scaling function itself since the integral in ( |5.11| ) depends just on 
The same holds, in fact, for all convolution operators 0]. 



Remark. The autocorrelation function of the scaling function (see ( |5.24|) ) has 2M — 1 
vanishing moments and its "zero moment" is equal to one (see ( p.25| ) and ( p.2(j| )). This 
implies that if we consider the representation of the derivative operator on the subspace 
Vq as a finite-difference scheme, such scheme has order 2M. For integral convolution 
operators, it implies that the asymptotics is accurate to order 2M (see ^ and below). 

The following proposition reduces the computation of the coefficients r/ to 
solving a system of linear algebraic equations. 

1. If the integrals in ( |5.10| ) or (|5.11|) exist, then the coefficients r^, / G Z in ( |5.10|) 
satisfy the following system of linear algebraic equations 



ri 



L/2 



^21 + 2 ^2k-l{r 2l-2k+l + f 2l+2k- 



2 ^ 

k=l 



and 



(5.12) 



(5.13) 



where 



L-2k 

0'2k-l = 2 ^ hihi+2k-l, 
i=0 



k = l,...,L/2 



(5.14) 



are the autocorrelation coefficients of the ffiter H. 

2. If M > 2, then equations (|5.12|) and ( |5.13D have a unique solution with a finite 
number of non-zero r/, namely, 7^ for —L + 2 < / < L — 2 and 



ri = — r_ 



(5.15) 



Solving equations ( |5.12| ), ( |5.13| ), we present the results for Daubechies' wavelets 
with M = 2, 3. For further examples we refer to [Q. 
1. M = 2 



and 



_ 9 _ 1 



2 1 

^^ = -3' ^^=12' 



We note, that the coefficients (—1/12, 2/3, 0, —2/3, 1/12) of this example can be found in 
many books on numerical analysis as a choice of coefficients for numerical differentiation. 
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Figure 6: Non-zero entries of a derivative operator in the non-standard form 



2. M = 3 

_ 75 _ 25 _ 3 

- 64' - ~128' - 128' 

and 

272 53 16 1 

ri = , To = , = , Ta = . 

365 365 1095 2920 

The structure of non-standard and standard forms of derivative operators is illus- 
trated in Figures |^ and 0. 

For the coefficients rj-"''' of d^/dx", n > 1, the system of linear algebraic equations 
is similar to that for the coefficients of d/dx. This system (and ( |5.12| )) may be written 
in terms of 

r(0=E^f"^e"«, (5.16) 



as 



f(0 = 2" ( |mo(e/2)|2f(e/2) + |mo(e/2 + vr)|2r(e/2 + vr)) , (5.17) 
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Figure 7: Non-zero entries of a derivative operator in the standard form 



where mo is the 27r-periodic function in ( |1.12|) . Considering the operator Mq on 27r- 
periodic functions 



(Mo/)(0 = |mo(e/2)|V(^/2) + |mo(e/2 + Tr)!^ /(^/2 + n), 



(5.18) 



we rewrite (|5.17| ) as 



Mof 



(5.19) 



so that f is an eigenvector of the operator Mq corresponding to the eigenvalue 2~". Thus, 
finding the representation of the derivatives in the wavelet basis is equivalent to finding 
trigonometric polynomial solutions of ( |5.19| ) and vice versa 

An important property of the wavelet representation of the (periodized) derivative 
operators (and, in general, pseudodifferential operators with homogeneous symbols) is 
that these operators have an explicit diagonal preconditioner in wavelet bases. 

We present here two tables illustrating such preconditioning applied to the stan- 
dard form of the second derivative. In the following examples the standard form of 
periodized second derivative D2 of size N x N, where iV = 2", is preconditioned in the 
wavelet basis by the diagonal matrix P, 



= PD^P, 



where Pn = 6ii2\ 1 < J < and where j is chosen depending on i,l so that 
A^/2^-^ + 1 < i, / < - A^/2-'', and Pnn = 2^". The matrix P is illustrated in Figure | 
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2 

2 

2 

2 

2 

2 

2 

2 

2 

2 
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4 

4 

4 

4 

4 

4 

4 

4 








8 

8 

8 

8 










16 

16 
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Figure 8: An example (n = 5) of the diagonal martix P used to rescale the matrix of 
the periodized second derivative in the wavelet system of coordinates. 

The following tables compare the original condition number k of and Kp of 



N 


K 




64 


0.14545E+04 


0.10792E+02 


128 


0.58181E+04 


0.11511E+02 


256 


0.23272E+05 


0.12091E+02 


512 


0.93089E+05 


0.12604E+02 


1024 


0.37236E+06 


0.13045E+02 



Table 2. 

Condition numbers of the matrix of periodized second derivative (with and without 
preconditioning) in the basis of Daubcchics' wavelets with three vanishing moments 

M = 3. 
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N 


K 




64 


0.10472E+04 


0.43542E+01 


128 


0.41886E+04 


0.43595E+01 


256 


0.16754E+05 


0.43620E+01 


512 


0.67018E+05 


0.43633E+01 


1024 


0.26807E+06 


0.43640E+01 



Table 3. 

Condition numbers of the matrix of periodized second derivative (with and without 
preconditioning) in the basis of Daubechies' wavelets with six vanishing moments 

M = 6. 



Fractional derivatives 

First, let us consider convolution operator T and the infinite matrix t^^_i ^\ i,l E 7i, 
representing Pj_iTPj_i on the subspace Vj_i. To compute the representation of PjTPj, 
we have (see e.g., formula (3.26) of [§) 

t?^ = j:j:h,h^i]-i^. (5.20) 

k=0 m=0 



It easily reduces to 



L/2 

=^21 + (^2k-l (t2i_2k+l + t2l+2k-l)- (5.21) 



fc=0 



where the coefficients 02^-1 are given in ( [5.141 ). 
We also have 



00 /■+00 



= / / K{x-y) d^dy^ (5-22) 

J— 00 J —00 

and, by changing the order of integration, we obtain 

tp) = 2W K{2^{l-y))<l>{y)dy, (5.23) 

where $ is the autocorrelation function of the scaling function ip, 

r+00 

$(?/) = / ip{x) ip{x - y) dx. (5.24) 
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It is easy to verify (see [§) that 



Hy)dy = 1, (5.25) 

J — oo 

and 

/+00 
<^{y)dy = 0, for 1 < m < 2M - 1. (5.26) 
-oo 

The vanishing moments of the autocorrelation function $ allow us to compute the el- 
ements of the matrix t^'' for large / and sufficiently fine scales j < 0. Expanding the 
kernel K in the Taylor series, we obtain from (|5.23|) 



t!^) = 2^K{2n) + ^^p^y^ y_ K^'''\2^{1 - y)) ^{y) dy, (5.27) 

where y = y{y, I) and TT^^m) denotes the (2M)th derivative ofK. The decay of /s:(2Af)(^2-'(/- 
y) ) for large / is faster than that of the original kernel (see ( [4.1| ) and ( |4.2| ) with an appro- 
priate choice of M) and ( [5.27| ) implies a one-point quadrature formula tp'' ~ 2^K{2H) 
for large / and sufficiently fine scales j < 0. 

Computing representations of convolution operators simplifies further if the sym- 
bol of the operator is homogeneous of some degree. Let us illustrate this using example 
of fractional derivatives. Defining fractional derivatives as 

^ f{y)dy, (5.28) 

where we consider 1,2.... IfQ;<0, then (|5.28|) defines fractional anti-derivatives. 
The representation of (9° on Vq is determined by the coefficients 

/+00 

if{x - I) (9» (x) dx, I e Z, (5.29) 
-00 

provided that this integral exists. 

The non-standard form = {Aj, Bj, Tj}j^z is computed via Aj = 2~°'^Aq, Bj = 
2~"-^i?o, and Tj = 2~°''To, where matrix elements and 7j_; of Aq, Bq, and Fq 
are obtained from the coefficients r/, 

L-l L-l 

ai = 2°'Yl 9kgk'r2i+k-k', (5.30) 

k=0 k'=0 

L-l L-l 

= 2" ^ gkhk'r2i+k~k', (5.31) 

fc=0 k'=0 
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and 



L-l L-1 

7i = 2" ^ J2 hk Qk' r2i+k-k' 

k=0 k'=0 



(5.32) 



It easy to verify that the coefficients r/ satisfy the following system of linear 
algebraic equations 



n = 2^ 



L/2 

^2/ + I X! ^2k-l{r 2l-2k+l + r 2l+2k-l] 
k=l 



(5.33) 



where the coefficients a2fc~i are given in (|5.14|) . Using (|5.271) , we obtain the asymptotics 
of r; for large I, 



ri 
ri 



1 1 1 



2M ' 







for / > 0, 

for / < 0. 



(5.34) 
(5.35) 



Example. 

We compute the coefficients r; of the fractional derivative with a = 0.5 for 
Daubechies' wavelets with six vanishing moments with accuracy 10~^. The coefficients 
for r;, / > 14 or / < — 7 are obtained using asymptotics 



ri 
ri 



for / > 0, 



for / < 0. 



(5.36) 
(5.37) 





Coefficients 




Coefficients 


I 


n 


I 


n 


-7 


-2.82831017E-06 


4 


-2.77955293E-02 


-6 


-1.68623867E-06 


5 


-2.61324170E-02 


-5 


4.45847796E-04 


6 


-1.91718816E-02 


-4 


-4.34633415E-03 


7 


-1.52272841E-02 


-3 


2.28821728E-02 


8 


-1.24667403E-02 


-2 


-8.49883759E-02 


9 


-1.04479500E-02 


-1 


0.27799963 


10 


-8.92061945E-03 





0.84681966 


11 


-7.73225246E-03 


1 


-0.69847577 


12 


-6.78614593E-03 


2 


2.36400139E-02 


13 


-6.01838599E-03 


3 


-8.97463780E-02 


14 


-5.38521459E-03 



Table 5. The coefficients {ri}i, I = —7, . . . , 14 of the fractional derivative a = 0.5 for 
Daubechies' wavelet with six vanishing moments. 
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VI Multiplication of matrices and fast iterative con- 
struction of the generalized inverse 

The standard and non-standard forms may be multiplied in fast manner if the matrices 
represent Calderon-Zygmund or pseudo-differential operators. Multiplication of matri- 
ces in the standard form is a straightforward algorithm |P, and requires at most 
0(A^log A^) operations. The algorithm for the multiplication of matrices in the non- 
standard form has been outlined in 0] and requires 0{N) operations. This is a signifi- 
cant improvement over 0{N^) operations for dense matrices which arise in the ordinary 
discretization of the operators from these classes. 

Fast multiplication algorithms give a second life to a great number of iterative 
algorithms. Indeed, powers of matrices maybe computed and so are other functions of 
matrices. Let us consider an iterative construction of the generalized inverse. In order 
to construct the generalized inverse A'^ of the matrix A, we use the following result ||31||: 



Let 0"! be the largest singular value of the mx n matrix A. Consider the sequence 
of matrices Xk 

Xk+i = 2Xfe — XkAXk (6.1) 

with 

Xo = aA\ (6.2) 

where A* is the adjoint matrix and a is chosen so that the largest eigenvalue of aA*A is 
less than one. Then the sequence Xk converges to the generalized inverse Al 



Combining this iteration with fast multiplication algorithms, we obtain an algo- 
rithm for constructing the generalized inverse in at most 0(iVlog^ iVlogi?) operations, 
where R is the condition number of the matrix. (By the condition number we understand 
the ratio of the largest singular value to the smallest singular value above the threshold 
of accuracy). 

The details of this algorithm (in the context of computing in wavelet bases) will be 
described in We note that throughout the iteration (|6.1D, it is necessary to maintain 
the "finger" band structure of the standard form of matrices Xk. Hence, the standard 
form of both the operator and its generalized inverse must admit such structure. We 
note that the pseudo-differential operators satisfy this condition. 

The following table contains timings and accuracy comparison of the construc- 
tionof the generalized inverse via the singular value decomposition (SVD), which is 0{N^) 
procedure, and via the iteration (|6.1| )- (|6.2| ) in the wavelet basis using Fast Wavelet Trans- 
form (FWT). The computations were performed on Sun Sparc workstation and we used 
a routine from LINPACK for computing the singular value decomposition. For tests we 
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used the following full rank matrix 



1 i = j 



where j = 1, . . . , iV. The accuracy theshold was set to 10 '^^ i.e., entries of Xk below 
10~^ were systematically removed after each iteration. 



Size N -xN 
128 X 128 
256 X 256 



SVD 
20.27 sec. 
144.43 sec. 



512 X 512 1,155 sec. (est.) 
1024 X 1024 9,244 sec. (est.) 



FWT Generalized Inverse L2-Error 
25.89 sec. 3.1 • 10"^ 

77.98 sec. 3.42 • lO"'* 



242.84 sec. 
657.09 sec. 



6.0-10-4 
7.7-10-4 



2^5 X 2^ 



9.6 years (est.) 



1 day (est.) 



We note that the iteration in (|6.1|) also allows us to compute the projector on the 
null space (see 0] for this and several other examples). 

The algorithm for the exponential is based on the identity 



exp(y4) = exp(2"-^y4) 



(6.3) 



First, exp(2-^y4) is computed by, for example, using the Taylor series. The number L is 
chosen so that the largest singular value of 2~^A is less than one. At the second stage 
of the algorithm the matrix 2~^A is squared L times to obtain the result. Similarly, sine 
and cosine of a matrix can be computed using the elementary double-angle formulas. 
Unlike the algorithm for the generalized inverse, this algorithm is not self- correcting. 
Thus, it is necessary to maintain sufficient accuracy initially so as to obtain the desired 
accuracy after all the mutiplications have been performed. 
Finally, as an example, let us consider the matrix 



/ -2 1 
1 -2 




1 





V 


















\ 




1 





-2 1 
1 -2 J 



(6.4) 
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Figure 9: Matrix Db computed via iterative algorithm of this Section with diagonal 
rescaling. Entries with the absolute value greater than 10~^ are shown black. 

which arises in the finite-difference formulation of the two-point boundary value problem. 
We note that the inverse of this matrix is sparse in the wavelet basis. As an illustration, 
in Figure ^ we display the inverse matrix Db~^ computed starting with matrix Db. 
Using the diagonal preconditioning (see Figure |), this computation involves only well- 
conditioned matrices |5|. 
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